# Figure_A06.R

# Part of the replication archive for 
#
#   Bullock, John G. 2020. "Education and Attitudes toward Redistribution in
#   the United States." British Journal of Political Science 50.


# This figure plots the effects of attendance-law strictness on educational 
# attainment. The estimates are stratified by -parents'- levels of educational 
# attainment.  The point is to show that the effects of attendance-law 
# strictness are concentrated among those whose parents didn't graduate from 
# high school.

library(Bullock, lib.loc = c(.libPaths(), 'packageLibrary'))       # for qw()
library(grid)
library(lattice)
library(lmtest)        # coeftest() helps with clustered standard errors
library(multiwayvcov)  # for clustered standard errors with lm objects
source("IV_setup.R")
dirOutput    <- 'float_output/'  
filenameStem <- 'Figure_A06'
PDF_title    <- 'Parental Education Moderates the Effects of Compulsory Attendance Laws' 
xJitter      <- c(-.025, .025)  # amount to offset results to left or right 



##############################################################################
# ESTIMATE MODEL 
##############################################################################
firstStageMod <- getFirstStage(IVModelsEnv$eqwlth.mod1)
parMod        <- update(
  firstStageMod, 
  educ ~ . - CA.fac + CA.fac*(parents.onlyOneHS_noBA + parents.HSnoBA + parents.onlyOneBA + parents.BA_orMore))
parMod.df     <- na.omit(GSS.df[, all.vars(parMod)]) 
parMod.lm1    <- lm(parMod, data = parMod.df, singular.ok = FALSE)


# CLUSTER THE STANDARD ERRORS
# make the state*year cluster variables
parModClusterVars <- droplevels(parMod.df$stateYoung:factor(parMod.df$yearYoungNorm))

# estimate the clustered vcov matrix
parModClusteredVcovMatrix <- cluster.vcov(
  model   = parMod.lm1,
  cluster = parModClusterVars)

# report estimates with clustered SEs
parMod.ests <- coeftest(parMod.lm1, vcov. = parModClusteredVcovMatrix)
parMod.ests <- parMod.ests[grepl('^CA', rownames(parMod.ests)), ]
rownames(parMod.ests) <- qw(
  "CAmod 
   CAstrict 
   CAmod:ParentsOneHS 
   CAstrict:ParentsOneHS
   CAmod:ParentsHS 
   CAstrict:ParentsHS
   CAmod:ParentsOneCollege 
   CAstrict:ParentsOneCollege
   CAmod:ParentsCollege 
   CAstrict:ParentsCollege")


# simplify vcov matrix for SE calculation
tmp <- grepl('^CA', rownames(parModClusteredVcovMatrix))
parModClusteredVcovMatrixSimple <- parModClusteredVcovMatrix[tmp, tmp] 
rownames(parModClusteredVcovMatrixSimple) <- rownames(parMod.ests)
colnames(parModClusteredVcovMatrixSimple) <- rownames(parMod.ests)



##############################################################################
# MAKE DATA FRAME OF DATA TO PLOT 
##############################################################################
dataToPlot   <- data.frame(
  parentType    = ordered(
    x      = rep(qw('parentsNoHS parentsOneHS parentsHS parentsOneCollege parentsCollege'), each = 2),
    levels =     qw('parentsNoHS parentsOneHS parentsHS parentsOneCollege parentsCollege')),
  # parentType    = ordered(
  #   x      = rep(qw('parentsNoHS parentsOneHS parentsOneCollege parentsCollege'), each = 2),
  #   levels = qw('parentsNoHS parentsOneHS parentsOneCollege parentsCollege')),
  CSL_type      = rep(c('CA in 8-10', 'CA geq 11')),
  coef          = c(
    parMod.ests['CAmod',    'Estimate'],
    parMod.ests['CAstrict', 'Estimate'],
    parMod.ests['CAmod',    'Estimate'] + parMod.ests['CAmod:ParentsOneHS',         'Estimate'],
    parMod.ests['CAstrict', 'Estimate'] + parMod.ests['CAstrict:ParentsOneHS',      'Estimate'],
    parMod.ests['CAmod',    'Estimate'] + parMod.ests['CAmod:ParentsHS',            'Estimate'],
    parMod.ests['CAstrict', 'Estimate'] + parMod.ests['CAstrict:ParentsHS',         'Estimate'],
    parMod.ests['CAmod',    'Estimate'] + parMod.ests['CAmod:ParentsOneCollege',    'Estimate'],
    parMod.ests['CAstrict', 'Estimate'] + parMod.ests['CAstrict:ParentsOneCollege', 'Estimate'],
    parMod.ests['CAmod',    'Estimate'] + parMod.ests['CAmod:ParentsCollege',       'Estimate'],
    parMod.ests['CAstrict', 'Estimate'] + parMod.ests['CAstrict:ParentsCollege',    'Estimate']),
  SE            = c(

    # SEs for children with 0 parents who finished high school
    sqrt(parModClusteredVcovMatrixSimple['CAmod', 'CAmod']),  
    sqrt(parModClusteredVcovMatrixSimple['CAstrict', 'CAstrict']),  
    
    # SEs for children with exactly one parent who finished high school, 
    # 0 parents who finished college
    sqrt(
      parModClusteredVcovMatrixSimple['CAmod', 'CAmod'] +
      parModClusteredVcovMatrixSimple['CAmod:ParentsOneHS', 'CAmod:ParentsOneHS'] +
      2*parModClusteredVcovMatrixSimple['CAmod', 'CAmod:ParentsOneHS']),
    sqrt(
      parModClusteredVcovMatrixSimple['CAstrict', 'CAstrict'] +
        parModClusteredVcovMatrixSimple['CAstrict:ParentsOneHS', 'CAstrict:ParentsOneHS'] +
        2*parModClusteredVcovMatrixSimple['CAstrict', 'CAstrict:ParentsOneHS']),

    # SEs for children with two parents who finished high school, 0 parents  
    # who finished college
    sqrt(
      parModClusteredVcovMatrixSimple['CAmod', 'CAmod'] +
      parModClusteredVcovMatrixSimple['CAmod:ParentsHS', 'CAmod:ParentsHS'] +
      2*parModClusteredVcovMatrixSimple['CAmod', 'CAmod:ParentsHS']),
    sqrt(
      parModClusteredVcovMatrixSimple['CAstrict', 'CAstrict'] +
        parModClusteredVcovMatrixSimple['CAstrict:ParentsHS', 'CAstrict:ParentsHS'] +
        2*parModClusteredVcovMatrixSimple['CAstrict', 'CAstrict:ParentsHS']),

    # SEs for children with exactly one parent who finished college
    sqrt(
      parModClusteredVcovMatrixSimple['CAmod', 'CAmod'] +
        parModClusteredVcovMatrixSimple['CAmod:ParentsOneCollege', 'CAmod:ParentsOneCollege'] +
        2*parModClusteredVcovMatrixSimple['CAmod', 'CAmod:ParentsOneCollege']),
    sqrt(
      parModClusteredVcovMatrixSimple['CAstrict', 'CAstrict'] +
        parModClusteredVcovMatrixSimple['CAstrict:ParentsOneCollege', 'CAstrict:ParentsOneCollege'] +
        2*parModClusteredVcovMatrixSimple['CAstrict', 'CAstrict:ParentsOneCollege']),
  
    # SEs for children with two parents who finished college
    sqrt(
      parModClusteredVcovMatrixSimple['CAmod', 'CAmod'] +
        parModClusteredVcovMatrixSimple['CAmod:ParentsCollege', 'CAmod:ParentsCollege'] +
        2*parModClusteredVcovMatrixSimple['CAmod', 'CAmod:ParentsCollege']),
    sqrt(
      parModClusteredVcovMatrixSimple['CAstrict', 'CAstrict'] +
        parModClusteredVcovMatrixSimple['CAstrict:ParentsCollege', 'CAstrict:ParentsCollege'] +
        2*parModClusteredVcovMatrixSimple['CAstrict', 'CAstrict:ParentsCollege'])),

  X = rep(1:(nrow(parMod.ests)/2), each = 2) + xJitter)
dataToPlot$upper <- dataToPlot$coef + 1.96*dataToPlot$SE
dataToPlot$lower <- dataToPlot$coef - 1.96*dataToPlot$SE



##############################################################################
# PRELIMINARIES FOR PRINTING THE FIGURE 
##############################################################################
postscriptBackground <- 'transparent'
plotLineColor        <- c('black', 'gray60')  # 'blue3' is better than 'blue2' on screen but too close to black in print.  Use 'grey37' if greyscale.
plotLineWidth        <- c(1.5, 1.5)
plotLineType         <- c('solid', 'solid')
panelBorderCol       <- 'black'
panelBorderWidth     <- .3

baseCexSize          <- 1.0
yLabelTextSize       <- .850 * baseCexSize   
xAxisTextSize        <- .670 * baseCexSize  # formerly .696
yAxisTextSize        <- xAxisTextSize
axisTickSize         <- .4
                     
PS_width             <- 7
PS_height            <- 7
panelHeight          <- list(2.500, "inches")
panelWidth           <- list(5.430, "inches")

N_parentEduc <- c(
  sum(parMod.df$parents.onlyOneHS_noBA),
  sum(parMod.df$parents.HSnoBA),
  sum(parMod.df$parents.onlyOneBA),
  sum(parMod.df$parents.BA_orMore))
N_parentEduc <- c(
  nrow(parMod.df) - sum(N_parentEduc),
  N_parentEduc)

xAxis <- list(
  draw   = TRUE, 
  limits = c(.75, nlevels(dataToPlot$parentType)+.25),
  labels = c(
    'neither parent\ngraduated from\nhigh school',
    'one parent\ngraduated from\nhigh school,\nneither from college',
    'both parents\ngraduated from\n high school,\nneither from college',
    'one parent\ngraduated\nfrom college',
    'both parents\ngraduated\nfrom college'),
  at     = dataToPlot$X[seq(1, nrow(dataToPlot), by=2)],
  tck    = c(axisTickSize, 0), 
  col    = panelBorderCol, 
  cex    = xAxisTextSize,
  alternating = 1, 
  relation    = 'same', 
  axs         = 'i')  # axs='i' means that there is no padding around the xlimits  
xAxis$labels <- paste0(
  xAxis$labels, 
  '\n(N = ', 
  prettyNum(N_parentEduc, big.mark = ','),
  ')')

yAxis <- list(
  draw        = TRUE,
  limits      = c(-.700, 2.10),  # -.615 1.95
  labels      = qw("0 .5 1 1.5"),
  at          = seq(0, 1.5, by = .5),
  rot         = 0,
  tck         = c(axisTickSize, 0),
  col         = panelBorderCol,
  cex         = yAxisTextSize,
  alternating = c(1, 2),
  relation    = 'free',       
  axs         = 'i')



##############################################################################
# PANEL FUNCTION
##############################################################################
myPanel <- function(...) {
  panel.xyplot(...)
  
  # TICK MARKS ON RIGHT-HAND SIDE OF PANEL
  panel.axis(
    side        = "right",
    at          = yAxis$at,
    draw.labels = FALSE,
    half        = FALSE,
    tck         = axisTickSize)
  
  # ADD VERTICAL CI LINES
  grid.segments(
    x0  = dataToPlot$X,
    x1  = dataToPlot$X,
    y0  = dataToPlot$lower,
    y1  = dataToPlot$upper,
    gp  = gpar(
      col = rev(plotLineColor),
      lwd = plotLineWidth,
      lineend = 'square'),
    default.units = "native")
  
  # ADD LABELS FOR EACH LINE WITHIN THE PANEL 
  trellis.par.set("clip", list(panel="off", strip="off"))
  grid.text(
    label = 'moderate\nattendance law',
    x     = unit( 4.275, "native"),
    y     = unit(-0.225, "native"),
    just  = c("left", "top"),
    # x     = unit(1.800, "native"),
    # y     = unit(0.195, "native"),
    # just  = c("right", "top"),
    gp    = gpar(
      cex        = baseCexSize*.6, 
      col        = plotLineColor[2],
      lineheight = .90))
  grid.text(
    label = 'strict\nattendance\nlaw',
    x     = unit(1.100, "native"),
    y     = unit(1.925, "native"),
    just  = c("left", "top"),
    gp    = gpar(
      cex        = baseCexSize*.6, 
      col        = plotLineColor[1],
      lineheight = .90))
}



##############################################################################
# PRINT THE PDF FILE
##############################################################################
pdf(
  file   = paste0(dirOutput, filenameStem, '.pdf'), 
  width  = PS_width, 
  height = PS_height, 
  paper  = "special", 
  title  = PDF_title,
  bg     = postscriptBackground)
trellis.par.set(
  axis.components = list(
    left = list(pad1 = .5)),   # distance between axis ticks and tick labels
  axis.line = list(
    alpha = 1, 
    col   = panelBorderCol,    # to eliminate lattice panel border, set col = "white"
    lty   = 1, 
    lwd   = panelBorderWidth),  
  layout.widths = list(
    ylab.axis.padding = .80),  # width between axis labels and ylab
  par.ylab.text = list(cex = yLabelTextSize),
  superpose.line = list(
    alpha = c(1,1), 
    col   = plotLineColor, 
    lty   = plotLineType, 
    lwd   = plotLineWidth))

parentalModerationPlot <- xyplot(
  coef ~ X, 
  groups   = CSL_type,
  type     = 'l',
  data     = dataToPlot, 
  panel    = myPanel,
  ylim     = panelYlim,
  strip    = FALSE,
  scales   = list(x = xAxis, y = yAxis),
  col      = plotLineColor,
  xlab     = '',
  ylab     = 'additional years of education',
)
print(
  x            = parentalModerationPlot, 
  panel.width  = panelWidth, 
  panel.height = panelHeight)
dev.off()

if (!(Sys.which('pdfcrop')=='' | Sys.which('perl')=='')  |  'pdfcrop' %in% installed.packages()[, 'Package']) {  # if "pdfcrop" is installed 
  system(
    paste(
      if (Sys.info()['sysname']=='Windows') paste(Sys.getenv('Comspec'), '/c ') else "", 
      'pdfcrop', 
      paste0(dirOutput, filenameStem, '.pdf'), 
      paste0(dirOutput, filenameStem, '_crop.pdf')))
  file.remove(paste0(dirOutput, filenameStem, '.pdf'))
  if (interactive() & Sys.info()['sysname']=='Windows') shell.exec(paste0(normalizePath(dirOutput), filenameStem, '_crop.pdf'))
}